#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#include <math.h>
#include <pthread.h>
#define NUMTHREAD 10

// <editor-fold defaultstate="collapsed" desc="variaveis globais">
double S = 1; // valor da ação

double E = 1; //preço de exercício da opção

double r = 1; //taxa de juros livre de risco (SELIC)

double sigma = 1; //volatilidade da ação

double T = 1.0; // tempo de validade da opção

double t;

int M = 2; //numero de iterações// 

double mean1, stddev1, confwidth1, confmm1, confmax1;

double * trials; 

pthread_t threads[NUMTHREAD];

double medvec[NUMTHREAD];

double desvioParcialVec[NUMTHREAD];

const unsigned int A = 1664525, C = 1013904223;
unsigned int x;

// </editor-fold>

// <editor-fold defaultstate="collapsed" desc="funcoes do metodo">

unsigned int next(unsigned int x) {
    return (x * A + C);
}

double max(double a, double b) {
    if (a >= b) return a;
    return b;
}

//mean sequencial
//double mean(double vec[], int m) {
//    int i;
//    double sum = 0;
//
//    for (i = 0; i < m; i++) {
//        sum += vec[i];
//    }
//    return sum / (double) m;
//}


//double stddev(double vec[], int N, double media) {
//
//    double sum = 0;
//    int i;
//
//    for (i = 0; i < N; i++) {
//        sum += (vec[i] - media) * (vec[i] - media);
//    }
//
//    return (sqrt((1.0 / (double) (N - 1)) * sum));
//}
void * desvioParalelo(void *arg){
    
    long offset;
    offset = (long)arg;
    int start;
    int end;
    int i;
    int passo = M/NUMTHREAD;
    start = offset * passo;
    end = (offset + 1) * passo;
    double soma = 0;
    double somaInterna = 0;
    
    
    for (i = start; i < end; i++) {
        
        somaInterna = trials[i] - mean1;
        somaInterna * somaInterna;
        soma += somaInterna;
        
    }
    desvioParcialVec[offset] = soma;
   
}

double confwidth(double stddev, int M) {
    return ( (1.96 * stddev) / (double) sqrt((double) M));
}

double confmm(double media, double confwidth) {
    return (double) (media - confwidth);
}

double confmax(double media, double confwidth) {
    return (media + confwidth);
}// </editor-fold>

// <editor-fold defaultstate="collapsed" desc="funcao blackScholes">

void * blackScholes(void *arg){
    
    long offset;
    offset = (long)arg;
    int start;
    int end;
    int i;
    int passo = M/NUMTHREAD;
    start = offset * passo;
    end = (offset + 1) * passo;
    double soma = 0;

    if(offset == NUMTHREAD - 1) 
        end--;
     
        for (i = start; i < end; i++) {
            
            x = next(x);
            
            double randomNumber = (double) x / (1LL<<32);

            t = S * exp( (r - ( (0.5) * sigma * sigma))  * T + sigma * sqrt(T) * randomNumber);

            trials[i] = exp( (-1) * r * T ) * max(t - E, 0);
            soma += trials[i];
            
        }
        medvec[offset] = soma/(double)passo;
        printf("\n medvec = %lf", medvec[offset]);
}

// </editor-fold>


int main(void) {
    
    x = time(NULL);
    double soma = 0;
    void *status;

    // <editor-fold defaultstate="collapsed" desc="inicializar variaveis e leitura da entrada">
    long i, j;
    
    scanf("%lf", &S);
    scanf("%lf", &E);
    scanf("%lf", &r);
    scanf("%lf", &sigma);
    scanf("%lf", &T);
    scanf("%d", &M);

    printf("Entradas: %lf %lf %lf %lf %lf %d\n", S, E, r, sigma, T, M);
    
    trials = (double*) malloc( M * sizeof(double) );
    
    // </editor-fold>
    
    // <editor-fold defaultstate="collapsed" desc="blackScholes">

    
    for (i = 0; i < NUMTHREAD; i++) {
        
        pthread_create(&threads[i], NULL, blackScholes, (void *)i );
    }

    for (i = 0; i < NUMTHREAD; i++) {
        
        pthread_join(threads[i], &status);
    }
    
    for(i = 0; i < NUMTHREAD; i++){
        soma += medvec[i];
        
    }
    
    mean1 = soma/NUMTHREAD;

    // <editor-fold defaultstate="collapsed" desc="calculo do desvio padrao">
    for (i = 0; i < NUMTHREAD; i++) {

        pthread_create(&threads[i], NULL, desvioParalelo, (void *) i);
    }

    for (i = 0; i < NUMTHREAD; i++) {

        pthread_join(threads[i], &status);
    }// </editor-fold>

    for(i = 0; i < NUMTHREAD; i++){
        soma += desvioParcialVec[i];
    }
    stddev1 = sqrt(soma/(double)(NUMTHREAD - 1));
    //stddev1 = stddev(trials, M, mean1);
    confwidth1 = confwidth(stddev1, M);
    confmm1 = confmm(mean1, confwidth1);
    confmax1 = confmax(mean1, confwidth1); 
    
    printf("S         %lf\n", S);
    printf("E         %lf\n", E);
    printf("r         %lf\n", r);
    printf("sigma     %lf\n", sigma);
    printf("T         %lf\n", T);
    printf("M         %d\n", M);
    printf("Confidence interval: (%lf, %lf)\n", confmm1, confmax1);
    
    // </editor-fold>

    free(trials);

    return (EXIT_SUCCESS);
}


